#include<stdio.h>
#include<cuda.h>
#include<time.h>

__global__ void MatrixMulKernel(float *Md,float *Nd,float *Pd,int Width)
{
	int tx=blockIdx.x*Width + threadIdx.x;
	int ty=blockIdx.y*Width + threadIdx.y;
	
	float Pvalue=0;
	
	for (int k=0;k<Width;++k) {
	
		float Mdelement=Md[ty*Width+k];
		float Ndelement=Nd[Width*k+tx];
		Pvalue += Mdelement*Ndelement;
	}
	
	Pd[ty*Width+tx]=Pvalue;
}

void HostMatrixMul(float *M,float *N,float *P,int Width)
{
	for (int i = 0; i < Width; ++i) {
		for (int j = 0; j < Width; ++j) {			double sum = 0;
			for (int k = 0; k < Width; ++k) {
				double a = M[i * Width + k];
				double b = N[k * Width + j];
				sum += a * b;
			}
			P[i * Width + j] = (float)sum;
		}
	}
}

int main(void)
{
	clock_t start=clock();
	cuInit(0);
	printf("\nInitialize CUDA %f", ((double)clock() - start) / CLOCKS_PER_SEC);
	
	int WIDTH=1000;
	size_t size=WIDTH*WIDTH*sizeof(float);
	int k;
	printf("\nInitialize Arrays %f", ((double)clock() - start) / CLOCKS_PER_SEC);
	
	start=clock();
	float *M = (float*)malloc(size);
	float *N = (float*)malloc(size);
	float *P = (float*)malloc(size);
	printf("\nHost Malloc %f", ((double)clock() - start) / CLOCKS_PER_SEC);
	
	start=clock();
	for(k=0;k<WIDTH*WIDTH;++k) {
		M[k] = rand() / (float)RAND_MAX;
		N[k] = rand() / (float)RAND_MAX;
		P[k] = rand() / (float)RAND_MAX;
	}
	printf("\nAssign Values %f", ((double)clock() - start) / CLOCKS_PER_SEC);
	
	float *Md, *Nd, *Pd;
	start=clock();
	cudaMalloc((void**)&Md, size);
	cudaMalloc((void**)&Nd, size);
	cudaMalloc((void**)&Pd, size);
	printf("\nCUDA malloc %f", ((double)clock() - start) / CLOCKS_PER_SEC);
	
	start=clock();
	cudaMemcpy(Md, M, size, cudaMemcpyHostToDevice);
	cudaMemcpy(Nd, N, size, cudaMemcpyHostToDevice);
	printf("\nMemcpy %f", ((double)clock() - start) / CLOCKS_PER_SEC);
	
	dim3 dimBlock(16,16);
	dim3 dimGrid(WIDTH/16,WIDTH/16);
	
	start=clock();
	MatrixMulKernel<<<dimGrid,dimBlock>>>(Md,Nd,Pd,WIDTH);
	printf("\nExecute Kernel %f", ((double)clock() - start) / CLOCKS_PER_SEC);
	
	start=clock();
	cudaMemcpy(P, Pd, size, cudaMemcpyDeviceToHost);
	printf("\nCopy Back %f", ((double)clock() - start) / CLOCKS_PER_SEC);
	
	start=clock();
	HostMatrixMul(M,N,P,Width);
	printf("\nExecute on Host %f", ((double)clock() - start) / CLOCKS_PER_SEC);
	
/*	printf("\n");*/
/*	for (i=0;i<WIDTH;++i)*/
/*	{*/
/*		for (j=0;j<WIDTH;++j)*/
/*			printf("%f ",P[i*WIDTH+j]);*/
/*		*/
/*		printf("\n");*/
/*	}*/
	
}
